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Abstract 



When the natural frequencies are allocated symmetrically in the Kuramoto model there exists an invariant torus of 
dimension [A'^/2] + 1 (A'^ is the population size). A global phase shift invariance allows to reduce the model to A'^ — 1 
dimensions using the phase differences, and doing so the invariant torus becomes [A'^/2]-dimensional. By means of 
perturbative calculations based on the renormalization group technique, we show that this torus is asymptotically 
stable at small coupling if is odd. If A'^ is even the torus can be stable or unstable depending on the natural 
frequencies, and both possibilities persist in the small coupling limit. 
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1. Introduction 

The Kuramoto model [l|2|3j has become the basic framework for the description of macroscopic synchro- 
nization; a phenomenon observed in a variety of natural and artificial systems [415j . Kuramoto [I] considered 
a population of all-to-all weakly coupled oscillators such that their interaction could be reduced to their 
phases: 



where 9j and ojj are, respectively, the phase and the natural frequency of the j-th oscillator, and e is the 
coupling strength. Kuramoto adopted a sinusoidal coupling function /(•) = sin(-) together with a symmetric 
frequency distribution of the natural frequencies, what resulted very useful for the theoretical analysis of 
the model. 

Originally, it was useful and instructive to consider the thermodynamic limit of ([1]), iV ^ oo. Finite-size 
effects have remained unsolved for a long time and only recently significative advances have been achieved 
^ 7 8 9 lOj. Also some attention has been recently devoted to the small- behavior of the Kuramoto model 
by Maistrenko and coworkers from the point of view of dynamical systems theory |11|12|13] . 
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In this paper we study the Kuramoto model with a finite population, and with the natural frequencies 
allocated symmetrically around the mean frequency. One of the reasons that motivates this problem is the 
fact that most works on the Kuramoto model have assumed that the natural frequencies are distributed 
according to a symmetric probability density, and as a consequence it is usual that numerical simulations 
are carried out selecting frequencies not at random, but reflecting the inherent symmetry of the frequency 
distribution. In particular, several works |ll|12|13|14j have recently investigated phase diagrams of the 
Kuramoto model with a finite population A'' under the assumption that the natural frequencies are allocated 
symmetrically. It has been shown that under these assumptions, finite N and symmetry of the natural 
frequencies, model ([T]) exhibits a peculiar type of chaos dubbed 'phase chaos'. 

Of more importance for this work is the finding in [Mj that the phase space contains an [A^/2]-dimensional 
invariant torus. This torus has been thought to be unstable {i.e. repelling) when e ^ [12|13|14] . This belief 
is probably motivated by the difficulty of investigating numerically the phase diagram for small e due to the 
extremely weak (in)stability of invariant sets in that limit. One of the purposes of this paper is to reveal the 
phase diagrams of the Kuramoto model at small coupling and different values of N. Our analytical results 
are obtained by using the renormalization group (RG) method, which is one of the singular perturbation 
methods. Our results firmly establish the stability of the mentioned invariant torus on rigorous mathematical 
grounds. In particular, we give general results for any finite N, and some small populations [N = 3, ..,7) 
are investigated in more detail. 

2. Basic definitions 

Due to the mean-field character of the Kuramoto model we may arbitrarily label the natural frequencies 
from the smaller to the larger: ivi < L02 < ■ • • < lo^. Moreover, by going into a suitable rotating framework 
we set the mean frequency equal to zero without loss of generality. For the discussion to follow it is worth 
to note that if there are not coincident natural frequencies {i.e. oji < L02 < ■ ■ ■ < wjv), some degree of 
synchronization, ^^i^ — ^^i^ ^"-"^ some (or all) i 7^ j, is only achieved for a coupling strength larger than 
some positive constant Sc. 

Next, we rewrite the Kuramoto model for the convenience of our analysis. Due to the invariance of the 
global quantity = X^il^i^'' ^^'^ Kuramoto model can be reduced in one dimension by changing to a 
new set of coordinates: Lpj — Oj+i — 9j, j = 1, ■ ■ ■ , N — 1. It is also useful to define "frequency gaps" 
Aj = i^j+i ^ j ~ ■ ■ ■ , N — 1. In the new variables the Kuramoto model has this structure 

<f3 ^ A-i +eEj{ipi,ip2,---,(PN-i), j = !,■■■, (2) 

Maistrenko and coworkers |12I13I14] found that if the natural frequencies are symmetrically selected 
{ui = —u)M-i+i Aj — Ajv_j) the Kuramoto model has an invariant manifold AA, namely a torus of 
dimension [N/2], which is given by = {(pi = (pjv-i}- (In the original coordinates, M is any of the tori 
that — parameterized by the invariant 8 — foliate an {[N/2] + l)-dimensional torus; thus if one selects = 0, 

M^{e,^ -eN-^+l}.) 

3. Main results 

In Refs. [12113114] it was shown that for some values of the parameters and e, the dynamics approaches 
the invariant torus M or an attractor A <Z M. They reported that the regions of parameters with stable M 
or A are close to, or inside, synchronization regions (what implies that e is larger than some positive number 
if uji ^ Thus, there is a common belief that far from synchronization, as the coupling strength goes 

to zero (e — > 0), the dynamics fills the whole phase space, the — 1 dimensional torus {^Pj}- So far, many 
numerical results confirmed this expectation [12113114] . However, we show in this paper that this is not true. 
We have found that in the limit e — > 0: 

(i) If N is odd the (iV — l)/2-dimensional invariant torus M is, unless there are special resonant conditions 
among the natural frequencies, asymptotically stable. 
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(ii) If N is even the iV/2-dimensional invariant torus Ai can be stable or unstable depending on the 
particular disposition of the natural frequencies. 

Statements (i) and (ii) are consequence of three Theorems to be proved by using the RG method, which is a 
powerful singular perturbation method for differential equations proposed in [15|16j . Recently, mathematical 
foundation of the RG method was given in |17|18j showing that the RG method is useful as well to investigate 
existence and stability of invariant manifolds. We present a brief review of the RG method in Section 7.1, 
and the proofs of the main theorems in this paper are shown in Section 7.2 and 7.3. 

4. Odd N 

In this section, we investigate the stability of the invariant torus M for odd N. For two particular cases, 
= 3 and 5, the phase diagrams are completely uncovered. 
One of the main theorems in this paper is as follows: 

Theorem 1. Suppose that N = 2M — 1 is an odd number. If the natural frequencies satisfy the following 
nonresonance condition: 

LOi = LOj if and only if i = j, 

LOk + LOj = 2uJi if and only if i = k = j or j = 2M — k,i ~ M, 
\ oji + ujj — ujk + oJi if and only if i = j = k = loij = 2AI — i,l = 2M — k, (3) 
3uji — bjj + Wfe + if and only if i ^ j = k ^ I, 

LOi + 2ujk — LOj + 2loi if and only if i = j, k = I or j = 2M ^ i,k = M, I ^ i, 

then there exists a positive constant eo: which depends on the natural frequencies, such that if < e < eo; 
the invariant torus A4 is asymptotically stable and the transverse Lyapunov exponents of A4 are of O(e^). 

Equation ([3]) can be rewritten as a condition for A^'s by using the relation uij = — X^^lj^ ^k- The proof 
of this theorem is given in Sec. 7.2. 

As £ — > 0, [A^/ 2] -frequency quasiperiodic dynamics on Ai is stable for almost all {A;}. Parameter regions 
on which ipi^s are (partially) phase- locked on A4 are very narrow. On such regions, there exist a fc-dimensional 
stable torus on Ai filled by fc-frequency quasiperiodic orbits (fc < [7V/2]). Further, we can prove that regions 
with phase- locking are narrower when the nonresonance condition is fulfilled than when it is not. 

Below we test the validity of Theorem 1 for iV = 3, 5, and 7. Focusing on particular cases will allow 
us to understand better how Theorem 1 applies in practical terms. For instance, for TV = 5 we make a 
complete analysis of the stability of M , showing what happens when the natural frequencies do not satisfy 
the nonresonance condition of Theorem 1. 



4.1. iV 3 



For A'^ = 3 it is particularly simple to prove the stability of Ai using basic theory of dynamical systems. 
The ODEs ruling the dynamics in {(pi} coordinates [Eq. ^] are: 



ifi = A + — [sin (^2 - 2sinv3i - sm{ipi + (^2)] 
1^2 = A + [sin tpi-2 sin tp2 - sin{ipi + (P2)] 



(4) 



where we are already assuming the symmetry Ai ~ A2 = A. The dynamics inside the invariant 1-torus Ai 
{i.e. a circle defined by tpi = ip2 = Lp) obeys 

if = A - -[sinip + sin{2ip)] (5) 
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and a transverse perturbation Sip = Lpi~ ip\ is governed by 

5iy9 = — ecos(/7(5(y9 + 0((Siy9^). (6) 
The transverse Lyapunov exponent (TLE) is hence: 

27r 

A_L = — e / P (if) cos (f dip, (7) 







with P{(p) — C/(p for £ smaller than the synchronization threshold Ec ~ 1.704A. C is a normalization 
constant such that ^ -P(v?) d(p — 1. Making an expansion of P{f) in terms of the small quantity e/A, it 
turns out that Aj_ becomes negative with a cubic dependence on e: 



1 ^3 

. 
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This result agreef0 with Theorem 1. Finally, it must be noted that the invariant torus exists for any odd 
interacting function f{(p) = —f{—ip) and not only for the particular choice f{(p) = sin((^). The transverse 
Lyapunov exponent is then: 
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Aj- - J /'M [2/(^)/(2¥') + f{2ip)] dip + 0{e'), (9) 



with J'{lp) EE df{(p)/dip. 
4.2. N = 5 

If TV = 5, Eq. (HD reads 

ipi = Ai + ^ 1^-2 sin (ySi + sin ip2 - sin(<y9i + 1P2) + sin((^2 + ^3) 

- sin(</7i + ip2 + ipa) + sin{ip2 + ips + ipi) - sin((^i + (p2 + ips + fiij , 
1^2 = A2 + (^-2 sin ip2 + sin ipi + sin (ps 

- sin((pi + (^2) - sin((/72 + fs) + sin((^3 + (P4) - sin((^2 + (fis + ^4)) , 

• A e / . . . ^ (10) 

1^3 = A2 + — (^-2 sm (y53 + sm (^4 + sm ip2 

- sin((p3 + (^4) - sin((/?2 + (^53) + sin((/?i + (/32) - sin((/3i + ip2 + ipsij , 
1^4 = Ai + -^(^-2sin(y54 + sin (^3 - sin(<y93 + ipi) + sin((^2 + (^3) 

- sin((/72 + </53 + </'4) + sin((^i + (^2 + '/'a) - sin((^i + (p2 + ips + ^iij ■ 

Again, symmetry is assumed (Ai = A4, A2 = A3), and hence there is an invariant 2-torus M — {ipi = 
f4, '■P2 — Va}- Since M is 2-dimensional, dynamics on M is nontrivial. Depending on values of Ai, A2 and 
e, the asymptotic dynamics on M. can be quasiperiodic or periodic (fixed points only exist above a finite e 
value unless Ai = A2 = 0). Quasiperiodic motion is generic when e — > 0. Periodic motion exist inside open 
sets in the phase diagram (so-called Arnold tongues) whose widths shrink to zero as e ^ 0. Inside an Arnold 
tongue there is (at least) one pair of stable-unstable (along the torus surface) periodic orbits whose average 
frequencies are related by a rational rotation number: {ipi) : {1P2) = n : m. We call a periodic orbit of that 



^ The nonresonance condition is not fulfilled if and only if A = 0, and in that case the TLE has a linear dependence on e: 
Ax = —s. Indeed, it is easy to see that the fixed point = of Eq. JSJ is stable for any e > 0. 
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(a) (b) (c) (d') (e'/f)(f) (g) (h) 
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Fig. 1. A schematic view of the phase diagram for N = 5 and small e. Inside the hatched regions A4 is unstable; and the dotted 
regions intend to represent finitely many disjoint hatched regions. Shaded regions are Arnold tongues with rotation numbers 
n:m = l:2, 1:1, 2:1 and 3 : 1 from left to right. 

type an n : TO locking solution. Arnold tongues touch the axis e = at the points where the ratio Ai : A2 is 
rational. Major Arnold tongues are born at (rational) Ai : A2 ratios corresponding to frequencies that do 
not fulfill the nonresonance condition of Theorem 1 . 

In what follows, we assumj^ A2 7^ 0. By rescaling time and e, we are allowed to divide Eq. PH)) by 
A2 to assume that A2 = 1 without loss of generality. Then, the nonresonance condition for = 5 gives 
Al 7^ 0, 1/2, 1, 2, 3, 4. We can prove the next theorem. 

Theorem 2. There exists a non-negative number £0 — £o(Ai) such that the invariant torus Ai is asymp- 
totically stable if < e < £o- eo(Ai) tends to zero as Ai 1/2, 1, 2, 3. 



A sketch of the strategy for proving this theorem is given in Sec. 7.3. And in Sec. 7.4 we show explicitly 
how the proof yields the phase diagram near Ai = 1/2 (the most intricate case). We note that Th. 2 asserts 
that the cases Ai = and 4 do not yield transversal instability of A4 despite of violating the nonresonance 
condition of Th. 1. A schematic view of the phase diagram of Eq. (fTO| for small e is represented in Fig. [U in 
which the invariant torus M is unstable in the tongue-shaped hatched regions. The n : to locking solutions 
for n:TO = l:2, 1:1, 2:1 and 3 : 1 exist in the gray regions. In particular, in the gray regions that arc not 
hatched there are transversally stable n : m locking solutions. Asymptotic expansions with respect to e/N of 
boundaries (a) to (h) and (a') to (h') in Fig.[l]are shown in the Appendix (the expansion for each line is done 
up to an order that completely unfolds the phase diagram). In two dotted regions emerging from Ai = 1/2, 
there are many disjoint unstable tongue-shaped regions, while exactly one unstable region emerges from each 
of the other resonances: Ai = 1, 2, 3. The existence of such many unstable regions emerging from Ai = 1/2 
is shown in Section 7.3 resorting to the RG method and with aid of numerical simulations. 

Next, we present numerical results corroborating our theoretical results for N — 5. The dynamics of 
infinitesimal perturbations transversal to the torus (Sipi = ipi — 1^94, 6ip2 = ip2~ ^3) is governed by two linear 
equations: 



^fi = -^cos(<^i 

Sip2 = --^C0S(<^1 



■ (p2){l + 4:COS(p2)S(pi + — [C0S(/J2 - COs{ipi 

ip2)sin^{tp2/2)S(fi ~ — [3cos</J2 + 2cos(</ji 



h ip2))] 5ip2- 



(11) 



From these equations we calculate the TLEs {\^^' > \^^') using the popular method by Benettin et al 



W>Af ) 

In Fig. Ufa) we plot the TLEs for two values of e observing that: (i) TLEs are almost everywhere negative 
as expected from Theorem 1, (ii) TLEs for e = 0.2 and 0.4 collapse when divided by e'^, except if one 
enters in a locking region nor fulfilling the nonresonance condition in Th. 1 (0 : 1 locking in the leftmost 



In the degenerate case A2 = the nonresonance condition of Th. 1 is violated if and only if uJi 



this case, Eq. JTDJ has a stable fixed point 1/3^ = 0, i = 1, . . . , 4 with TLEs A 



(1) 



A 



(2) 



: (Al = A2 = 0). In 
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Fig. 2. (a) Transverse Lyapunov exponents scaled by as a function of Ai for e = 0.2,0.4 (A'^ = 5, A2 = 1). (b-e) Regions 
where the largest TLE becomes positive. (It is instructive to compare the panel (b) with the result in Fig. obtained by 
the RG method.) 
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Fig. 3. TLEs as a function of e for N = 5, A2 = 1 and Ai = (1 + \/5)/2 (the golden mean). 

part of the panel). Figure [3] shows a log-log plot of the TLEs as a function of e for a specific value of Ai 
(arbitrarily chosen with the constraint that the nonresonance condition is satisfied). We find a nice power 
law A\_' ^ = -|AY' 'I oc e'^, as expected from Theorem 1. In Figs. [2l(b-e) we depict the largest TLE in four 
regions about lockings 1:2, 1:1, 2:1, and 3:1, finding that it becomes positive in short intervals as 
advanced in Theorem 2. These intervals match with the analytical expressions in Appendix. 



4.3. = 7 



For N = 7 the system's dimension is too large to perform a detailed analysis around all relevant resonances. 
But still the dimension is small enough to carry out intensive numerical simulations. We fixed A3 = 1 and 
measured the largest TLE at different values of Ai and A2. The initial condition in Ai was random (what 
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3A^=A3 

0.5 1 1.5 , 2 2.5 3 

^1 

Fig. 4. A1-A2 plane for N = 7 with A3 = 1 and e = 0.3{Af/2) min; A^. In the green regions the dynamics in Ai is transversally 
stable (i.e. there is an attractor A. C Ai). Inside the black regions Ai is transversally unstable. (Red) lines indicate the loci of 
the most important resonances among {A^}. 

should not be a problem if as we expect multistability is not common at small e). In our simulations the 
coupling strength was e = 0.3(A^/2) min^ A.;. This value of e is small enough to ensure the systems is far 
from synchronization, and at the same time large enough to make convergence times not exceedingly long for 
our computational resources. We may expect this value of e to capture the fundamental phenomenology as 
e ^ 0. The result is presented in Fig.[51 and shows that the dynamics on Ai is transversally stable in almost 
all the A1-A2 plane, except close to some resonances. These resonances should correspond to combinations 
of Ai's not fulfilling the nonresonance condition of Th. 1 (but also one may not exclude finite-e effects). 
The result is very much equivalent to the result for N — 5, but with unstable regions organized around 
resonances involving three instead of two frequencies. 



5. Even TV 

In the symmetric Kuramoto model at small coupling the invariant torus M is almost always stable if the 
population size is odd. However, we show in this section that if N is an even number M can be both stable 
and unstable in large regions of the parameter space spanned by the natural frequencies. We analyze the 
case iV = 4 in detail by means of the renormalization group (RG) method, and the case iV = 6 is studied 
using numerical calculations. Both cases share features that should be common to any even number A'^ > 4. 



5.1. = 4 

When = 4, Eq. (O with Ai = A3 is written as 
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G) (k') (1) (m') (n) (o')^p,^ (q) (r') (s) 



(t) 



(i') (i) G')^ I (k) (l'))|(m)(n')^ (p)) f(q') (r) \ |(s')^^\ 



(u) 




Fig. 5. A schematic view of the phase diagram of Eq. I|12|l for small e. There are no attractors A C JVl in the hatched regions, 
and gray regions indicate the Arnold tongues. 



Lpi = Ai + — 2 sin (ySi + sin Lp2 

- sin(<y9i + (^2) + sin(<y92 + V's) - sin((^i + <^2 + V'a)) , 
(/'2 = A2 + — 2sin(y52 + sin(^i + sin (^3 — sin((^i + (^2) — sin((^2 + (/^a) 



(12) 



"^3 = Ai + — ( -2 sin + sin (^2 



■ sin((/72 + (^3) + sin((/7i + (/32) - sin((^i + (/32 + V3)) ■ 



In this case, the invariant torus is given by the 2-dimcnsional torus Ai = {<pi = (/33,(p2}- Note that there 
wiU exist n : m lockings on M Uke for N — 5. 

The Unear ODE governing infinitesimal deviations off the torus {5ip = ipi ~ (p^) is: 



Sip = -- [c0S(pi + C0s((pi + ip2)] Sip, 



(13) 



and it determines the TLE. 

In what follows, we suppos^ A2 7^ 0. By dividing Eq. (fT2|) by A2, we can assume A2 — 1 without loss 
of generality. 



Theorem 3. There exists a non-negative number Eq = eo(Ai) such that if < e < eo, the invariant 
torus Ai is asymptotically stable for Ai > 1 and unstable for Ai < 1. eo(Ai) tends to zero as Ai 
0, 1/3, 1/2, 1, 2, 3, 4. The transverse Lyapunov exponent of A4 is of O(e^) if < e < Eq- 

A sketch of the strategy for proving this theorem is given in Sec. 7.3, while the detailed calculation is 
omitted. A schematic view of the phase diagram of Eq. for small e is depicted in Fig. O in which there 
are no attractors ^ C in the hatched regions, and the n : m locking solutions for n:m = 0:l,l:3, 1: 
2, 1 : 1, 2 : 1, 3 : 1 and 4 : 1 exist in the gray regions. In particular in the regions that are gray but not 
hatched, there are stable periodic orbits on A4. Asymptotic expansions with respect to e/N of boundary 
curves (i) to (u) and (i') to (u') in Fig. [5] are shown in Appendix. 

We numerically calculated the TLE from Eq. ([T3|), and Figs. [6] and [7] demonstrate that the TLE is 0{e^) 
as stated in Theorem 3. Ai is mainly stable when Ai > 1. Like for = 5 there are "switching tongues" in 
which the stability of A4 is different from the dominant stability in its neighborhood. We do not show them 

^ If A2 = there are two situations: (i) If Ai = 0, Eq. I|12|l has a stable fixed point tpi = 0, {i = 1, 2, 3) with = — e. (ii) 
If Ai ^0 there is a stable periodic orbit (1:0 locking solution). Dynamics for ip2 is well described by averaging the second 
equation of Eq. I|12| l with respect to ipi and 933: 932 = — 2-^sin</32- It proves that 1^2 = is stable. Substituting (/P2 = into 
Eq. 1 II2I I. we obtain the equation for ip = {pi = ifj,: tp = Ai — (£/A')(2 sin + sin(2(p)). The transverse Lyapunov exponent is 
calculated using Eq. Mil in the same way as that of Af = 3, obtaining = -£^/(16A2). 
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Fig. 6. Transverse Lyapunov exponent scaled by e as a function of Ai for two values of e {N = 4, A2 = 1). The scale of the 
y-axis is nonlinear to better discern the sign of the TLB. There is a perfect overlap of both data sets (e = 0.2 and 0.4), except 
at frequency lockings : 1 and 1 : 1. The stability change at Ai = 1+ is very sharp but continuous. 
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Fig. 7. TLEs as a function of e for Af = 4, A2 = 1 and (a) Ai = (-l+V5)/2 (the inverse of the golden mean) (b) Ai = (l + V5)/2 
(the golden mean). For the case (a) the invariant torus }A is unstable while it is stable for (b) as stated in Theorem 3. 

in Fig. [6] not to overwhelm the reader with details. We note nevertheless that the analytical expressions (see 
Appendix) of boundary curves (i) to (u) have been corroborated by our numerical simulations. 



5.2. TV > 6 



We resort to extensive numerical simulations to study the case iV = 6, as we did already in Sec. 14.31 for 
N = 1. Figure [8] summarizes the result of our simulations. There are regions with stable A^, and regions 
with unstable Ai . Some low order resonances give rise to the border among these regions (analogously to the 
Ai = A2 border in the = 4 case). Other resonances give rise to thin strips where the stability switches 
(again in good analogy to the = 4 case). 

If Ai is sufficiently large, oscillators 61 and 6q rotate so fast (|wi| ^ I1JJ2.3I) that their influence on the other 
oscillators averages out. In this regime, (in) stability of the invariant torus is ruled by the four oscillators 
with the central frequencies. It may be perceived in Fig. [5] that for Ai 3> A2.3 there is a switch of stability 
at A2 ^ A3 = 1 in qualitative agreement with the Ai = A2 +0(e^) border in the A^ = 4 case (the transition 
is not closer to A2 = 1 due to the finiteness of e and Ai). 

This simple argument based on the averaging method can be extended to an arbitrary (even) population 
size: In the limit Ai 00, Eq. ^ for A^ = 2M + 2 becomes equivalent to Eq. ([2|) for A^ = 2M because 
influence of ipi and (p2A/+i on the other oscillators averages out. Thus by induction, we conclude that the 
situation of A^ = 4 is typical for general even N; that is, regions of frequency space where A4 is stable and 
regions where it is unstable coexist in parameter space, none of them disappearing as £ ^ 0. This property 
quite differs from the phase diagrams for odd A^. 
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Fig. 8. A1-A2 plane for N = 6 with A3 = 1 and e = 0.3{N/2) min^ A^. In the green regions the dynamics in Ai is transversally 
stable (i.e. there is an attractor A. C J\4). Inside the black regions Ai is unstable. (Red) lines indicate the loci of the most 
important resonances among {A^}. In the regions with small Ai or A2 the convergence was too slow to guarantee accurate 
results (we have indications that typically the TLBs are 0(e'') as for A'' = 4). 

6. Discussion 

One of the important results of our paper is the rather surprising fact that quahtative properties of the 
Kuramoto model (with symmetrically allocated natural frequencies) depend crucially on whether N is odd 
or even. We establish precise mathematical criteria for the stability of the invariant torus . It is remarkable 
that in many cases this torus is asymptotically stable at arbitrarily small coupling. For = 4, 5, we have 
completely uncovered stability changes caused by resonances. 

The renormalization group method has been successfully used in this paper. In the literature, first- and 
second-order RG equations have been employed for constructing approximate solutions to weakly perturbed 
ODEs. In this paper, we have used this technique in the Kuramoto model up to third and fifth order for 
odd N and for A'^ = 4, respectively. RG equations are quite helpful for studying the stability of invariant 
manifolds, and they provide as well the orders of magnitude of the transverse Lyapunov exponents. 

In the context of coupled oscillators, the phenomenon known as 'phase chaos' consist in the appearance 
of a high-dimensional chaotic attractor [20121114] due to the interaction of phase variables (neutrally stable 
variables in the uncoupled limit). For the Kuramoto model it has been reported in '14] that (with a symmetric 
allocation of the natural frequencies) when increasing the coupling from zero the Lyapunov exponents split 
from a degenerate set at zero to a set of [N/2] — 1 positive, [N/2] — 1 negative, and 2 (3 if odd) zero 
Lyapunov exponents. In this paper we show that the invariant torus M is often stable in the small coupling 
limit. In contrast it is not proven yet that phase chaos indeed persists in the £ — s- limit. Here we note that 
our Theorem 1 applies in a finite range < e < Eq irrespective of how large N is. Notice nevertheless that 
Eo — > as ^ 00. 

Other implications of our paper refer to numerical simulations. In this context our paper is particularly 
relevant because the Kuramoto model has been usually considered together with a symmetric frequency 
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distribution of the natural frequencies, and in turn simulations, with a finite population size, are often 
carried out selecting frequencies that reflect the inherent symmetry of the frequency density. Our results 
also evidence the important role of resonances among the natural frequencies for the stability of the invariant 
torus A4. This must serve as a warning about the risks of using highly resonant frequencies when tackling 
generic properties of the system. In particular, it has become popular to consider evenly spaced natural 
frequencies (A^ = Aj for all which is probably the most resonant case. 

We may conjecture that the results reported here for the Kuramoto model (with a sinusoidal coupling 
function) should also be observed in a family of odd coupling functions. In fact, in the iV = 3 case, see 
Eq. ([9]), a family of functions shares stability and scaling of the transverse Lyapunov exponent. Nonetheless 
in what concerns the unstable regions (tongue-shaped for N = b) and their scaling one may expect important 
differences depending on the coupling function /. (We suspect this should be the case because of the similarity 
with the phase-locking regions, which depend on the harmonics of the interaction function [22j.) 

Finally, note that under a small enough symmetry-breaking perturbation A4 will get deformed into an 
invariant torus A4' with the same stability. Therefore our results may apply to situations where the symmetry 
is weakly broken. 

7. Outline of the proofs of theorems 
7.1. Brief review of the RG method 

The renormalization group (RG) method is one of the singular perturbation methods for differential equa- 
tions which provides approximate solutions as well as approximate invariant manifolds and their stability. 
Recently, it is shown that the RG method unifies and extends traditional singular perturbation methods, such 
as the averaging method, the multi-time scale method, the normal forms theory and so on. In this section, 
we give a brief review of the RG method following [17118] to prove Theorems 1 in the next subsection. 

Consider the system of differential equations on a compact manifold M of the form 



where e G M is a small parameter. For this system, we make the following assumption (A): 

(A) The vector fields gi{t,x), i = 1,2, - • • are with respect to time t G R and C°° with respect to 

X € M. Further, gi are almost periodic functions with respect to t uniformly in x € M , the set of whose 

Fourier exponents has no accumulation points on M. 

In the case of the Kuramoto model H]), M is an iV-dimensional torus. Note that under the change of 

coordinates 9j = Xj + ujjt and ipj = Xj + Ajt systems ([1]) and ([2]), respectively, are transformed into the 

form of Eq. (fT4|l with gi = for i > 2, and satisfying the assumption (A). 

Substitute x = xq + exi + 8^x2 + • ■ ■ into the right hand side of Eq. \1A\ and expand it with respect to e. 

We write the resultant as 



^x = egi{t,x) +e'^g2{t,x) +e^g3{t,x) H , x e M, 



(14) 



00 



00 



s'^gkit, xo + exi + 8^X2 + 



) ^ e''Gk{t, xo,xi, - ■ ■ , Xk-i). 



(15) 



k=l 



For instance, Gi, G2 and G3 are given by 



G2{t,xo,xi) 



Gi{t,xo) 




(16) 
(17) 



Gsit, Xo,Xi,X2) 



(18) 



respectively. With these G^'s, we define the C°° maps : M M to be 
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i?i(y)= lim - fGiis,y)ds, (19) 

t 

u^l\y)= f{Gi{s,y)-Ri{y))ds, (20) 



and 



R,{y) = ^lim - {G,{s,y,ui'Hy), ■ ■ ■ ^u^t'Hy)) ^ -^{y)R^-k{y))ds, (21) 

k=l y 

i i-1 (fc) 

u'j^\y) = J (g,(s, u^^\y), • • • , ^.(^-^^(y)) - - i?,(2/))ds, (22) 



for i = 2, 3, • • •, respectively, where J* denotes the indefinite integral, whose integral constants are fixed 
arbitrarily. We can prove that Ri are well-defined {i.e. the limits exist) and u^^^ are bounded in i e M. Along 
with Ri and u[^\ we define the m-th order RG equation for Eq. p4p to be 

y = eRiiy) + e^R2iy) + ■■■ + e'^Rmiy), y & M, (23) 

and the m-th order RG transformation al™"* to be 

at\y) = y + eu^^\y) + ■■■+ e"^u["'\y), y e M. (24) 

Roughly speaking, we can show that the m-th order RG transformation x — a["^^ (y) brings the system 
((T4| into the system of the form y — eRi{y) + • • • + e"''Rm{y) + e™^^5(t, x, e), where S is bounded in t e K. 
It means that the m-th order RG equation is e™+^-close to the original system (fT4|) and thus it is useful to 
construct the flow of (|14p approximately. Since the RG equation is an autonomous system while the original 
system ([T^ is not, to analyze the RG equation is easier than that of the original system. The next theorem 
is one of the fundamental theorems of the RG method. 

Theorem A [17,18j. Suppose that Ri{y) ~ ■ ■ ■ ~ Rk-i{y) = and £^Rk{y) is the first non-zero term 
in the RG equation. If the vector field Rk{y) has a boundaryless compact normally hyperbolic invariant 
manifold A/", then for sufficiently small e > 0, Eq. (fT4|) has an invariant manifold M^, which is diffcomorphic 
to N . In particular, stability of A^ coincides with that of A/". 

This theorem is used to investigate the stability of the invariant torus M. and the n : m locking solutions 
of the Kuramoto model. 



7.2. Proof of Theorem 1 

In this section, we give the proof of Theorem 1. 

Proof of Theorem 1. Suppose that N = 2M — 1 is an odd number and w.^'s are allocated symmetrically 
as was assumed. Put 6i = Xi + LOit and rewrite Eq. ^ of the form of Eq. (fTi]) . If the natural frequencies 
satisfy the nonresonance condition ([3]), its third-order RG equation is given by 
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VM 



^ \ sin(2yM - Vk - y2M-k), 



/ 1 1 

= 1^ 1^1^ ^°^^yi 



22/M + 2/2M-i) 



E 



\k^i,2M-i 



-2 E 

k=£i,2M-i 

+2 E 

- E 



1 



1 



■ sin(?;i - j/fe - y2M-k + y2M-%) 



s\n{yi -yk- y2M-k + y2M-i) 



(25) 



LUk[i^i - i^ki 
1 



k^MA,2M-i 



sin(2j/A/ - 2/fc - y2M-k) 

sin(?/j - yfc - 2/2Af-fc + y2M-i) 



-2 



E 



k^M,2M-i 



sin(?;i - 2?;m + y2M-i) , (« 7^ M)- 



Note that the first order term vanishes and the expansion begins with the second order term. Since the 
invariant torus M corresponds to the solution + y2M-i = c (constant), we put (j)i ~ yi + j/2M-i and 
4>M — '2yM- Then we obtain the system of (jji : 



64£^ 



9a/ 



32e-^ 



M-1 ^ 

E-^ sin((/)M - ^fe), 
/, I. 

A/-1 



iV3 ^ W 

fc=l 



1 

2 ' 

M-1 



hi) - 4 J2 



k^i 



uJt - ojZ 



■ sm 



(26) 



+4E 



3 sin((/)M - 0fc) , (i = 1, • • • , M - 1). 



k^i « K / 

Now that the second order term vanishes and Theorem A for fc = 3 is apphcable to this system. We can 
prove that the eigenvalues of the Jacobian matrix of the r.h.s. of ([26l) at the fixed point 0^ = c (i = 1, • • • , M) 
have negative real parts, except a zero eigenvalue that results from the rotation invariance of Eq. ^ (or the 
degree of freedom of the constant c). A proof of this fact is outlined as follows: 

Let be the Jacobian matrix of the r.h.s. of ([^5]) at the fixed point = c {i — 1, ■ ■ ■ , M). By using the 
cofactor expansion, it is easy to show that the characteristic polynomial of is calculated as 



det (A/ 



iV3 

32^ 



J) = A-det (A/ + ^Af-i), 



(27) 



where the matrix -4a/ -i is given as 
/ 3 ^ 4 



Am- 



1 ~ ^2 
A/-1 



a;5 



2 3,-^ 

1 1 2 j,_^2 



^1 ~ "^M-l 

4 



— Lot 



-'A/-1 



^A/-l 



OJ 



A/-1 



M-1 

E 



■'M-l 



-'A/-1 



(28) 



M-l fc^A/-l '^M-l "^k j 
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Eq. (|27p shows that J has a zero eigenvalue A = 0. Now it is sufhcient to prove that aU eigenvalues di Am-\ 
have positive real parts. To prove it, let 

be the characteristic equation AeXiXl — Am-\) = of Am-\- We show the inequalities ■ ■ ■ < 

and /ifli'^\ /l/li^^ • ■ ■ > by induction on AI. Since f^-i invariant under the permutation of 
ai, • • • , QfAf-i, we can show that /V-i the form 

Jm-1 - °M-1 Z2~2 72 ' "M-1 e ^- I'JlJj 

Since Am-i \ | as aM-i ^ 0, ^ /if -2^ as um-i 0. Now induction on M proves the 

desired inequalities. 

Thus, the solution (pi = yi + y2M-i = c (i = 1, • • • , M) of the RG equation is asymptotically stable and 
this proves that the invariant torus M is asymptotically stable for small £ > 0. Note that the degree of 
freedom of c does not appear in the Lpj coordinates [Eq. ■ 

If = 3(Af — 2), the nonresonance condition ([3]) is reduced to Ai 7^ and Theorem 1 recovers the results 
obtained in Section 4.1. 




7.3. Sketch of the proofs of Theorems 2 and 3 



Theorems 2 and 3 are also proved by using the RG method though we need much harder analysis to 
obtain asymptotic expansions of the boundary lines in Figs. [T] and [S] In this section, we offer the strategy 
to prove Theorems and to derive the asymptotic expansions, which is also valid for any N . 

In what follows, we assume A2 = 1 in Eqs. (fTO|) . (fT2|) as was mentioned. Our strategy to prove Theorems 
2 and 3, and to obtain boundary lines, is summarized as follows: 

(i) Put (pi^ Xi + Ait and rewrite Eqs. and into the form of Eq. (fT4)) . 

(ii) Derive the RG equations up to third-order for = 5 and to fifth-order for N = 4. The forms of RG 
equations depend on Ai . Find the set of values Ai , which gives the nonresonance condition, at which 
RG equations take different forms from the others. We find that the nonresonance conditions are given 
by Ai ^ 0,1/2,1,2,3,4 for TV = 5 and Ai ^ 0, 1/3, 1/2, 1, 2, 3, 4 for = 4. 

(iii) Investigate the stability of the invariant torus for the RG equation satisfying the nonresonance condition 
as was done in Sec. 7.2. In this step, the proof of Theorems 2 and 3 ends. 

(iv) To find the Arnold tongues in Figs. [1] and El let cq be a resonance value obtained in step (ii) . Put 

Ai=co + cie/N + C2e^/N^ + ■■■ (31) 

in Eqs. (fTO|) and (|12p. and derive the RG equations. 

(v) Investigate the stability of the invariant torus for the resultant RG equations and find values ci, C2, • • • 
at which the stability changes. Then, Eq. ([31]) gives an asymptotic expansion of a boundary line of the 
Arnold tongue emerging from Ai = cq in the phase diagram. 

In the next subsection, we calculate asymptotic expansions of boundary lines (a), (b), (a') and (b') in 
Fig. [T]to confirm this strategy for A'' = 5 and Ai near 1/2. Other expansions of boundary lines (c) to (u) 
and (c') to (u') in Figs. [1] and [5] are obtained analogously, and their derivation is omitted; the results are 
given in the Appendix. 



7.4. Phase diagram near Ai = 1/2 for N ^ 5 



In this section, we derive asymptotic expansions of lines (a) and (b) , as well as asymptotic expansions of 
lines (a') and (b'), which are boundaries of the 1 : 2 Arnold tongue for A^ = 5. We also show that there are 
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many disjoint unstable regions of A4 emerging from A1/A2 = 1/2 as is shown in Fig.[TJ 

To investigate the phase diagram of iV = 5 near A1/A2 = 1/2, put A2 = 1, Ai = l/2 + cie/N + C2e'^/N'^ + 
■ ■ ■ and put (fii = xi + 1/2, (p2 = X2 + 1, ip^ = X3 + 1 and (p4 ~ X4 + t/2 in Eq. (fTO|) . Then Eq. ((TO)l takes the 
form of Eq. (fT4|) and the RG method is appHcable. The third-order RG equation for the system is given by 











d 


2/2 


e 





di 


2/3 


^ N 














(-1 

2 



+ C2 - - cos(y2 - 2/3) + 2 cos(yi + V2 - Vs. - Vi) 
^ cos(?/2 - 2/3) 



V 



1 



10 
1 1 

10 ^ 2 



cos(2/2 - V-i) 



-7: + C2 - - cos(y2 - 2/3) + ^ cos(j/i + y2 - 2/3 - 2/4] 



£3 



(32) 



/ 26ci 2ci , N 1 . N \ 

~3~ + ^ cos(yi + 2/2 - 2/3 - 2/4) - 2 ^^"(^yi - j/2) » 

27 1 77 

sin(y2 - 2/3) - Y2 ^^^^^^2 ~ '^y^'' ~ 45 ^^'^'^^i + 2/2 - 2/3 - Vi) 

96Cl „ . /„ N 1 • /r, \ 1 • / \ 

- 2sm(2?/i - y2) + sm(2yi - 2/3) + ^ sm(y2 - 2/3) 

1 . / ^ N 3 . , „ N 26 . , 
+ T sm y2 - 2j/4) - ;r sm(y3 - 2^4) + — sm(?/i + 2/2 - 2/3 - 2/4] 
6 2 45 

96ci 3 . , N 1 • \ 1 • / \ 

'^5- + 2 sin(22/i - 2/2) - g sm(22/i " 2/3) - 2q sin(2/2 - 2/3) 

1 , , , , 26 , 

- sm(?/2 - 22/4) + 2 sm(2/3 - 22/4) - ^ sm(2A + 2/2 - 2/3 - 2/4) 

26ci 2ci , ^ 1 . 

~3~ ~ ~9~ '^"^^^i + 2/2 - 2/3 - 2/4) + sm(22/i - 2/3) 
27 1 77 

V +^sin(2;2 - 2/3) + 2^111(2/3 - 22/4) + Y^sin(2/i +2/2 - 2/3 - 2/4) ) 

This system has the solution 2/1 — Hi, 2/2 = 2/3: which corresponds to the invariant torus A4. The linearized 
equation for Eq. ([32]) along the solution 2/1 = 2/4 j 2/2 = 2/3 is given as 



dt 





3 / 




£3 


\sy2j 





II, , 551 5 cosy 

/__(44 + 15cosy) + ^ 

52 13 cos 2/ 113 15 cos 2/ 




(33) 



2 90 4 

where ^2/1 — Vi ~ 2/4: '^2/2 = 2/2 ^ 2/3: ^'^d where 2/ = 2/(^) is a solution of the equation 



dy 2e 



dt N 



119 



-rr = + 1:77 2c2 - — + 
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£3 /4664 



1 



7V3 V 225"^^ + '^^-4^^"^'' 



(34) 



which is obtained by putting 2yi — 2/2 = 2/i 2/4 = 2/i: 2/3 = 2/2 in Eq. ([32| . The existence of 1 : 2 locking 
corresponds to the existence of a stable fixed point of y, and the stability of the invariant torus A4 is maps 
to the stability of the trivial solution Syi — Sy2 = in Eq. ([55]) . Stability depends on the coefficients c,; in 
Eq. HMD. 

(i) When ci ^ 0, then we can apply the averaging method to Eqs. p3l [34)) . Averaged with respect to 2/, 
Eq. ([55]) is rewritten as 
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4 6 



- 0.05 




Fig. 9. (a) Lyapunov exponents (Ai > A2) of system 1331 1 forced by 1361 1 {e/N = 1 is arbitrarily adopted). In the interval 
0.125552... < C3 < 5.3695..., Ai is mostly positive except at hundreds of narrow windows where it becomes negative (to 
coincide with A2). (b) In this panel we see how Ai presents dips that accumulate at C3 = 1/8; the dips are so narrow that the 
sampling is not able to resolve the intervals where Ai = A2. 




7V3 



-154/45 -551/90 \ / 5yi 
52/45 113/90 j I (52/2 



(35) 



It is easy to verify that the trivial solution of this system is stable because the eigenvalues of the matrix 
in the right hand side are given by (— 13 ± \/231z)/12. It proves that the invariant torus M is stable 
if ci 7^ 0. 

(ii) When ci = and C2 ^ 119/30, we can apply the averaging again to obtain Eq. ([35)1 what proves that 
M is stable in the same way as (i). 

(iii) When ci = and C2 = 119/30, Eq. ^ becomes 



dt 



2e^ 

iV3 



C3 



— siny 



(36) 



(iii-a) If jcsl < 1/8, the above equation has a stable fixed point y — y* such that siny* — 8C3, cosy* = 
(1— 64c§)^/^. It corresponds to the 1 : 2 locking solution because y = 2yi— y2. The disappearance of 
the fixed point at jcaj ~ 1/8 marks the boundaries of the Arnold tongue: asymptotic expansions 
(a') and (b') in Appendix. It is easy to investigate the stability of the trivial solution of the 
linearized equation ()33|) with constant coefficients. Indeed, we can show that the trivial solution 
is unstable if and only if 



, , 1 / 8681209 \ ^ 

C3<- , = 0.124838. 

4 V 17391218 + 5894^8738809/ 



(37) 



And this proves that the invariant torus A4 is unstable in the region surrounded by the lines (a) 
and (b) given in Appendix, 
(iii-b) If jcal > 1/8, the linearized equation p3|) is a linear system with a time periodic coefficient. It 
is well known that stability of a trivial solution of such a system is determined by the Floquet 
exponents although wc can not calculate them analytically in general. We examine the stability 
of the trivial solution of Eq. (|33p by calculating the Lyapunov exponents numerically. 

Figure [5] shows that there are many disjoint intervals of C3 on which the trivial solution of 
Eq. (|33p is unstable. It proves that there are many disjoint unstable regions of the invariant 
torus A4 emerging from Ai = 1/2. These unstable regions are inside the region limited by Ai = 
i + + 0.125552;^ + 0{e^) and Ai = i + + 5.3695^ + 0{e^); and inside a twin 

region with opposite signs in the cubic terms (dotted regions in Fig.[T]). 
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Appendix : Asymptotic expansions of boundary curves in Fig. [T] (TV — 5) and Fig. [5] {N — 4) 
Asymptotic expansions in e of lines (a) to (h) in Fig. [1] which divide stable regions and unstable regions 
of the invariant torus Al for iV = 5, arc 

, ^ 1 119 £2 If 8681209 V'^ £3 4^ 

(a,b)Ai = - + — -^Tt L^...... . .... ....... J T7^ + 0{e% 



2 30 iV2 4 V 17391218 + 5894V8738809/ 
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. N A 1 2 £2 689009 £4 

= ^+3]V^- 362880 iV^ + ^^^ )' 

^ 3 ^ V 241920 V ^4 ^ V ^' 

. N A o 5 £2 139753 £'^ 

^^-2--^-^^ + 0(£^), 

(f) A.^2-^4+f-H^ + fc.)4 + 0(£^), 
^ ^ 3 7V2^V 48384 7 7V4 ^ ^ 

161 £2 / 191\/T848T- 24281 e 



3 



(^'^)^^-^-W]v^^l, 28800 ) iV^ + ^(^')' 

where fci and ^2 £^re some positive constant which are not obtained analytically. 

Asymptotic expansions of boundaries (a') to (h') of the n : m Arnold tongues in Fig. [T]are 

. A 1 2 £2 684751 ^, 

) = 1+ 3 iV^ - 241920 A^ + ^^^ )' 
/ ,/N A 1 5 £2 1027631 ^, 

^^+3iV^ + ^4192(7 A^ + ^(^^)' 
, „ , „ 5 £2 291373 £4 

^^-2-3iV^- 241920A^ + ^^^)' 
2 £2 136655 £4 
(f) ^i-2-3]^-l^]^ + 0(£ ), 

Asymptotic expansions of the lines (i) to (u) in Fig. [5l which divide stable regions and unstable regions 
of the invariant torus A4 for iV = 4, are 

+ 0(£2), 



(i) 


Ai 


(j) 


Ai 


(k) 


Ai 


(1) 


Ai 


(m) 


Ai 


(n) 


Ai 


(P) 


Ai 


(q) 


Ai 


(r) 


Ai 




Ai 



V805 iV 

1 193 £2 6767219 e'^ 102957740948201549 



3 30 84000 59535676200000 

I 193 £2 145573959475337 5 
3^"30"iV2" 3308633832000 
1 23 £2 2 £3 202777 72568 £^ g 

2^"6"iV2^ m ^ 7560 TV" ^ 1575 TV^ + 

1 593 £2 3^ 

2 + 106iV^ + ^(^)' 

4^2^ _ 24289 e* 5 

3 7V2 4320 
„ 168449 £2 
^-l750^iV^ + ^(^)' 

2 47 £2 1 £3 241399 1339 ^ 

^ 15 iV2 ^ 4 iV3 ^ 378000 ^ 56700 TV^ 

239 £2 168549449269 £4 5 
~ l2"iV2 ~ 41374295424 TV^ 

239 £2 11279063 £4 12019873522646189 s'^ 7 
^ l2"iV2 ^ 3259872 TV^ ~ 6354207699840000 TV^ ^ ' 
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, , ^ .719 1658470531 V63584702081 ^, 

(t, u) Ai = 4 T ^ + Oie^). 

^ ' ^ 90 7V2 265356000 TV^ ^ 3302208 ^ ' 

Asymptotic expansions of boundaries (i') to (u') of the n : m Arnold tongues in Fig. [S] are 



(0 


A 

Ai 




+ 0(6% 


(J ) 


A 

Ai 


1 

= 3^ 


193 


8249819 5 
84000 N-^ ^ ^ 


(k) 


A 

Ai 


1 

= 3^ 


193 


6767219 £" 1286745722182601 
84000 TV" 742041300000 


(1) 


A 

Ai 


1 

= 2^ 


23 £2 

6 
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